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1.  INTRODUCTION 


1 . 1  Problem  Overview 


The  capability  of  carrying  out  numerical  computations  of  flow  fields 

around  finned  missile  configurations  is  of  significant  importance  to  the  U.S. 

1 1-2  1 

Army.  '  Spatial  marching  techniques  have  previously  been  applied  to 

predict  three-dimensional  viscous  flow  about  axisymmetr ic  bodies  and  finned 

(  3-7  1 

geometries  of  interest  to  the  Army  for  low  angles  of  attack.  '  For  such  a 
technique  the  solution  is  marched  spatially  down  the  body  in  the  predominant 
flow  direction.  However,  this  technique  is  only  applicable  to  problems  for 
which  the  flow  is  supersonic  and  does  not  contain  imbedded  subsonic  regions  or 
regions  where  the  flow  separates  in  the  marching  direction.  A  time  marching 
subdomain  adaptive  technique  that  utilizes  Navier-Stokes  (NS)  equations  for 
separated  regions  or  in  the  vicinity  of  solid  surfaces,  and  Euler  equations 
elsewhere,  was  proposed  to  overcome  the  shortcomings  of  spatial  marching 

/  g  \ 

techniques.  '  This  technique  was  titled  Strategy  Combinina  P»gionally  Adapted 
Processes  (SCRAP). 

The  overall  stated  objective  of  the  research  is  to  carry  out  3-D 
numerical  simulations  of  finned  missiles  and  guided  projectiles  which  have  a 
gap  between  the  inner  fin  edge  and  the  body  surface.  The  simulation  of  the 
flow  field  about  these  configurations  at  a  nori-zeiu  angle  oi  dLtack  is 
important  from  the  standpoint  of  understanding  the  physical  phenomena  and 
predicting  the  magnitude  of  the  forces  acting  on  the  fins  and  the  projectiles. 
The  flow  field  simulation  of  the  full  configuration  (rather  than  just  the  fins 
and  gap  regions)  is  necessary  because  the  flow  field  associated  with  the  body 
affects  the  one  associated  with  the  fins  and  vice-versa  at  high  angles  of 
attack.  The  flow  fields  are  geometrically  complicated  because  of  the  3-D 
configurations  and  the  close  proximity  of  fins  and  body.  This  makes  the 
process  of  initial  grid  generation  complicated.  Also,  the  fin/body 
configurations  have  complex  physics  associated  with  them.  This  may  include 
shock  waves,  boundary  layers,  recirculation  regions,  etc.,  and  their 
interactions.  It  is  important  to  consider  numerical  schemes  that  have  good 
shock-capturing  properties.  Full  NS  computation  with  an  appropriate 


turbulence  model  may  be  prohibitively  expensive.  However,  full  NS  solutions 
is  only  needed  in  the  gap  regions  and  in  the  boundary  layers  or  recirculation 
regions.  In  addition,  it  may  also  be  needed  in  the  regions  of  interactions  of 
features  like  shock  waves  and  boundary  layers.  For  much  of  the  computational 
domain,  Euler  solution  is  adeguate.  Therefore,  it  makes  sense  to  consider 
different  numerical  schemes  in  different  selected  regions  to  speed  up  the 
convergence  of  the  calculation. 

Although  the  SCRAP  technique  was  originally  proposed  to  simulate 
transonic  or  supersonic  flow  about  finned  missiles  and  guided  projectiles, 
this  methodology  can  be  applied  to  any  physically  and/or  geometrically  complex 
flow  field.  The  issue  of  the  complexity  of  the  grid  generation  can  be 
partially  simplified  by  considering  a  multiblock  gridding  technique.  In  this 
approach  each  component  of  the  overall  configuration  has  an  independent  grid 
structure  which  is  called  a  block.  The  multiblock  SCRAP  technique  simplifies 
physical  complexity  by  associating  a  numerical  scheme  to  each  of  the  blocks. 
Thus,  one  can  have  viscous  blocks  where  NS  equations  are  solved  and  inviscid 
blocks  where  one  carries  out  the  solution  to  Euler  equations.  Such  a 
scheme-splitting  technique  can  efficiently  model  complex  physics  if  the 
domains  of  viscous  and  inviscid  regions  are  known  a  priori.  The  SCRAP 
technique  addresses  this  issue  by  considering  interaction  blocks  where  a 
decision  is  locally  made  whether  to  use  Euler  or  NS  equations.  This  dynamic 
physics  switching  is  termed  here  as  physics  (or  equation)  adaptation .  The 
interaction  blocks  will  be  important  for  the  cases  where  the  viscous/inviscid 
interaction  region  develops  as  part  of  the  numerical  solution.  The  SCRAP 
technique  also  implements  subdomain  adaptation  in  which  the  algorithm 
automatically  implements  additional  spatial  resolution  by  locally  dividing 
cells  in  the  regions  of  enhanced  gradients  of  selected  flow  field  variables. 

1 • 2  Overview  of  Phase  I  Accomplishments 

This  subsection  summarizes  the  key  accompl i shment s  of  Phase  I.  These 
accomplishments  pertain  to  the  implementation  of  the  overall  SCRAP  technique 
in  two  spatial  dimensions.  Extension  to  3-D  will  be  carried  out  in  Phase  II. 
The  accomplishments  include: 
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Development  of  an  Interactive  Block  Grid  Generator  ( I BGG ) . 


•  Development  of  a  general  purpose  code  that  implements  physics 
and  subdomain  adaptation  on  multiblock  grids. 

•  Demonstration  of  the  potential  of  SCRAP  technique  for  a  model 
problem. 

tor  the  Interactive  Block  Grid  Generator  ( I BGG )  the  flow  domain  is 
divided  into  blocks,  and  each  of  these  blocks  is  associated  with  a  given 
numerical  scheme.  These  blocks  are  classified  as  Euler,  Navier-Stokes , 
Thin-Layer  Navier-Stokes  (TLNS),  or  interaction,  depending  upon  the  scheme 
being  utilized.  Interaction  blocks  are  the  ones  in  which  a  local  decision  is 
made  for  each  of  its  cells  whether  to  use  an  Euler  or  NS  solver.  In  other 
words  physics  adaptation  is  carried  out  on  these  interaction  blocks.  The  user 
specifies  the  nodes  of  the  boundaries  of  the  blocks  which  must  have 
node-to-node  matching  at  inter-block  interfaces.  Utilizing  the  interblock 
boundary  point  information  the  interior  grid  of  each  of  the  blocks  is 
generated.  Subsequent  to  that,  all  the  nodes  that  are  multiply  defined  are 
removed  and  the  overall  boundary  nodes  are  adjusted.  Then,  all  the  cells  with 
the  same  numerical  algorithm  are  stored  together  in  the  computer  memory;  these 
cells  need  not  be  physically  contiguous.  Cells  in  each  block  can  be 
subsequently  spatially  divided  by  the  algorithm  to  incorporate  resolution 
adaptation.  Regions  of  cells  can  also  be  initially  manually  subdivided  to 
incorporate  additional  resolution  or  pre-embedding.  Various  stages  of  the 
grid  assembly  can  be  graphically  viewed  and  debugged  interactively  in  this 
I  BGG.  For  example,  nodes  can  be  moved  individually  or  in  concert  with  an 
elliptic  grid  generator.  Similarly  the  information  about  the  connectivity  of 
various  elements  (cells,  nodes,  boundary  points)  can  be  queried  interactively. 

The  overall  algorithm  is  decoupled  into  a  viscous  and  inviscid 
contribution.  This  means  that  Euler  equations  are  solved  everywhere, 
including  the  viscous  region.  This  is  followed  by  the  determination  of 
viscous  cells  in  the  interaction  blocks.  These  cells  are  lumped  together  with 
cells  in  the  viscous  regions.  Additional  viscous  contributions  (corresponding 
to  NS  or  TLNS  equations)  are  computed  for  these  viscous  cells  and  added  to  the 
inviscid  contributions.  The  decoupled  nature  of  the  algorithm  allows 
different  time  steps  to  be  taken  for  viscous  and  inviscid  contributions.  This 
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is  in  line  with  the  concept  of  local  time  stepping  for  carrying  out  solucions 
to  steady  state  problems.  Physics  adaptation  is  incorporated  by  finding  the 
maximum  shear  stress  in  the  cells  of  viscous  and  interaction  blocks.  Those 
cells  of  the  interaction  block  for  which  the  shear  stress  is  more  1%  of  the 
maximum  shear  stress  are  regarded  as  viscous  cells.  Subdomain  adaptation  is 
carried  out  by  examining  the  gradients  of  selected  flow  variables  in  the 
domains.  The  cells  for  which  the  first  differences  of  these  variables  exceed 
a  threshold  limit  are  locally  subdivided. 

The  potential  of  the  SCRAP  technique  is  demonstrated  for  a  supersonic 
flow  over  an  8%  circular  arc  cascade.  This  demonstration  clearly  establishes 
that  simultaneous  physics  and  resolution  adaptation  works  for  the  sample 
problem.  The  technique  provides  significant  savings  in  CPU  time  and  storage 
over  conventional  structured  grid  algorithms  employing  a  full  NS  solver.  An 
increase  of  about  two  orders  of  magnitude  in  computational  speed  was  achieved 
when  compared  to  a  conventional  full  NS  solution  over  a  fine  grid  of  the  same 
mesh  size  as  the  smallest  grid  size  chosen  by  the  spatial  resolution 
adaptation. 

1 . 3  Overview  of  Report 


The  following  section  summerizes  the  technical  objectives  of  Phase  I  and 
the  status  of  accomplishments.  Section  3  describes  the  details  of  the  methods 
utilized  to  achieve  each  part  of  the  technical  objectives.  Section  4 
describes  the  actual  numerical  results  obtained  by  utilizing  the  SCRAP 
technique.  Conclusions  of  this  report  are  summarized  in  Section  5. 
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2.  PHASE  I  TECHNICAL  OBJECTIVES 


The  stated  specific  objectives  of  the  Phase  I  research®  were: 

1.  Implement  integration  schemes  for  solving  Euler,  thin-layer 
Navier-Stokes  (TLNS)  and  complete  Navier-Stokes  (CNS)  equations 
within  different  relevant  zones. 

2.  Generate  data  management  connectivities  that  will  allow  viscous 
or  inviscid  formulation  on  an  element-by-element  basis. 

3.  Specify  conditions  to  be  applied  at  zonal  interfaces. 

4.  Demonstrate  the  potential  of  using  SCRAP  technique  for 
predicting  resultant  forces  on  a  practical  finned  missile 
configuration  with  a  gap  between  fin  and  body. 

5.  Document  the  results  and  conclusions  of  Phase  I  in  a  Final 
Technical  Report. 

All  these  Phase  I  objectives  have  been  successfully  completed  for  two  spatial 
dimensions.  Extension  to  three  spatial  dimensions  will  be  carried  out  in 
Phase  II.  Although  the  development  of  an  interactive  block  grid  generator 
( IBGG )  was  not  mentioned  in  the  list  of  objectives,  it  was  carried  out  as  part 
of  the  effort  to  facilitate  initial  grid  generation.  Similarly  the  end 
product  of  Phase  I  was  promised  to  be  a  feasibility  study  report  that  would 
establish  the  effectiveness  of  a  numerical  procedure  that  incorporates  physics 
adaptation  for  distinct  zones  of  a  multiblock  system.  We  have  established 
this  for  subdomain  adaptation  as  well.  The  methods  used  to  achieve  each  part 
of  the  Phase  I  technical  objectives  are  detailed  in  the  following  sections. 

The  Phase  I  effort  successfully  demonstrated  the  basic  feasibility  of 
simultaneous  physics  and  resolution  adaptation  for  multiblock  grids.  The 
technique  provides  significant  savings  in  CPU  time  and  storage  over 
conventional  structured  grid  algorithms  employing  a  full  NS  solver  over  a  fine 
grid  of  the  same  mesh  size  as  the  smallest  grid  size  chosen  by  the  spatial 
resolution  adaptation.  The  SCRAP  technique  yields  essentially  the  same 
results  as  the  global  calculation  over  a  fine  grid  system.  The  results  were 
also  compared  with  previous  computations  and  were  found  to  be  in  good 
agreement.  This  successful  demonstration  of  proof  of  concept  in  Phase  I  has 


5 


formed  the  basis  for  extension  to  3-D  in  Phase  II.  The  ultimate  objective  of 
the  proposed  effort  is  to  deliver  a  complete  and  general  purpose  SCRAP 
software  package  that  would  be  user  friendly,  modular,  and  would  incorporate 
all  the  relevant  physics  pertaining  to  complex  3-D  objects  in  flow  fields 
encompassing  subsonic  to  supersonic  Mach  numbers. 
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3.  NUMERICAL  DETAILS 


This  section  describes  the  details  of  the  methods  utilized  to  achieve 
each  part  of  the  Phase  I  technical  objectives.  The  section  starts  with  the 
description  of  a  general  set  of  full  Navier-Stokes  equations  which  can  be 
utilized  as  a  physical  model  for  the  flow  simulation.  Subsection  3.2 
describes  the  normalization  of  these  equations.  This  is  followed  by  the 
description  of  the  inviscid  and  viscous  terms  in  the  numerical  solver  in  terms 
of  an  unstructured  finite  volume  cell-vertex  scheme.  The  algorithm  is  shown 
to  be  decoupled  into  contributions  from  viscous  and  inviscid  terms. 
Subsection  3.4  describes  the  implementation  of  boundary  conditions. 
Subsection  3.5  presents  a  von  Neumann  stability  analysis  for  a  one-dimensional 
convection-diffusion  model  problem  which  suggests  that  different  time  steps 
can  be  used  for  the  contributions  from  viscous  and  inviscid  terms  in  the 
overall  algorithm.  Subsection  3.6  describes  an  interactive  block  grid 
generator  that  was  developed  as  part  of  the  Phase  I  effort.  Subsection  3.7 
describes  the  implementation  of  physics  adaptation  in  the  cells  of  interaction 
blocks. 

A  brief  description  of  the  new  contributions  and  their  significance  is  as 
follows: 

•  Development  of  an  Interactive  Block  Grid  Generator  that  allows 
grid  splitting  in  user-defined  blocks  corresponding  to  specific 
schemes.  The  interior  domain  block-interfaces  can  be  moved  to 
generate  "better"  grids. 

•  Decoupling  of  the  integration  scheme  into  contributions  from 
inviscid  and  viscous  terms.  This  allows  efficient 
implementation  of  the  algorithm  for  viscous  and  inviscid  blocks. 
Different  time  steps  for  these  terms  accelerate  the  convergence 
to  steady  state. 

•  Efficient  implementation  of  physics  adaptation  for  cells  of 
interaction  blocks.  For  these  cells  it  is  locally  decided 
whether  to  use  Euler  or  NS  solver. 

•  Development  of  a  general  purpose  code  that  efficiently 
implements  physics  and  subdomain  adaptation  on  multiblock  grids. 
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The  complete  Navier  Stokes  equations  describe  a  general  3-D  compressible 
viscous  flow  of  a  gas  and  are  given  by  the  following  vector  form 

2>U  SF1  SiG1  iiH1  2>FV  J>GV  2>HV 

-  +  -  +  -  +  -  =  -  +  -  +  -  ( 1 ) 

at  ax  ay  az  ax  ay  az  '  ' 

where  t  represents  the  time  and  (x,y,z)  the  spatial  coordinates  as  independent 
variables.  The  superscripts  i  and  v  signify  inviscid  and  viscous  terms, 
respectively.  The  state  vector  U,  inviscid  flux  vector  F1  and  viscous  flux 
vector  Fv  are  given  by 


u 

[P, 

pu,  pv,  pw,  e] 

T 

(2) 

F1  = 

[pu, 

,  pu2  +  p,  puv, 

puw ,  (e+p)u)T 

(3) 

Fv  = 

(0, 

Txx'  Txy '  Txz' 

UTxx+VTxy+WTxz-<*xlT 

(4) 

where  p,  p,  (u,v,w),  and  €  are  the  density,  preBBure,  velocity  components,  and 
total  energy  per  unit  volume,  respectively  and  T  denotes  the  transpose 
operation.  Furthermore  t  and  q  represent  the  stress  tensor  and  heat  flux, 
respectively.  The  representation  of  the  flux  vectors  G  and  H  is  analogous  to 
that  of  F  and  is  rot  shown  here.  Euler  equations  are  obtained  by  setting  the 
viscous  flux  vectors  equal  to  zero. 

Additional  relations  are  required  to  complete  the  system  of  equations. 
An  equation  of  state  that  connects  pressure  with  total  energy  per  unit  volume 
is 


£ 

P 


V l 
2 


+  C_.T  -  £ 

P  p 


E 

P 


(5) 


where  V  is  the  magnitude  of  total  velocity 


Cp  is  the  specific  heat  at  constant  pressure  and  T  is  the  temperature  of  the 
gas.  A  calorically  perfect  gas  is  assumed  for  which  the  specific  heat  Cp  is  a 
constant . 

The  shear  stress  tensor  t  is  usually  represented  as  a  linear 
phenomenological  law  in  terms  of  first  and  second  viscosity  coefficients  (u 
and  X) .  This  in  indicial  notation  is  represented  by 


T 


i-j 


M( 


+ 


+ 


S 


k=l 


auk 

J>xk 


(7) 


where  S^j  is  the  Kronecker  delta 


S^j  =1  if  i=j  ;  0  otherwise  .  (8) 

The  first  and  second  viscosity  coefficients  are  usually  connected  by 

utilizing  Stokes  hypothesis  which  states  that  the  sum  of  the  diagonal  elements 
of  the  shear  stress  tensor  must  equal  zero.  This  hypothesis  is  obviously 
valid  for  a  static  fluid  but  is  assumed  to  hold  even  for  a  fluid  in  motion  and 
results  in 

.  2 

X  =  -yy  .  (9) 

As  an  example,  the  viscous  shear  stresses  for  a  2-D  flow  are  given  by 


= 

2y ) 
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+  X( 

bv 

bw 

2 

bu 

bv 

T 

XX 

+ 

bx 

-  + 

^y 

bz  ) 

=  jy  (2 

bx 

*Y) 

(10) 

=  (X- 

2y ) 

bv 

+  X( 

bu 

bw 

2 

bv 

bu 

'T 

yy 

+ 

by 

—  + 
bx 

bz  * 

"  F(2 

by 

bx  * 

(ii) 

=  M  ( 

bu 

bv 

Txy 

by 

+ 

bx  * 

• 

(12) 

The  dynamic  viscosity  coefficient  y  is  a  strong  function  of  temperature  T  and 
is  frequently  modeled  by  the  Sutherland  law 


M 

Hr 


Tr  +  S 
^  T  +  S  ^ 


(13) 
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where  the  subscript  r  denotes  a  reference  condition  and  S  is  a  constant  which 
is  characteristic  of  a  given  gas.  For  air  at  moderate  temperatures,  S  =  110.4 
K. 

The  heat  flux  q  is  generally  related  to  the  thermal  conductivity  k  by  the 
Fourier  law 


q  =  -k  VT  =  -k  (—  i  +  —  j  +  —  k)  (14) 

A  constant  Prandtl  number  is  generally  utilized  as  a  relation  between  dynamic 
viscosity  and  thermal  conductivity 


For  air  at  normal  conditions,  Pr  =  0.72. 


3 . 2  Normalization  of  Governing  Equations 


Let  the  subscript  r  indicate  some  reference  conditions  and  denote  the 
non-dimensional  quantities  by  asterisks,  i.e.,  define 


X 

=  x*Lr 

y 

ii 

K 

* 

C 

z 

=  z*Lr 

t 

=  t*tr 

p 

* 

=  P  Pr 

G 

=  e*er 

u 

* 

=  u  Ur 

V 

* 

=  V  vr 

W 

* 

=  W  Wr 

p 

* 

=  P  Pr 

T 

=  T  Tr 

M 

★ 

=  M  Mr 

In  order  to  keep  the  form  of  the  dimensional  and  normalized  equations 
invariant,  the  continuity  equation  dictates 


(17) 


The  Inviscid  part  of  the  momentum  flux  terms  yield 


u 


2 

r 


£r 

Pr 


(18) 
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whereas  the  inviscid  part  of  energy  flux  term  implies 


=  Pr 


prur 


(19) 


The  viscous  flux  terms  of  the  momentum  equations  imply  the  following 
non-dimensional  parameter  of  the  fluid  flows 


Re 


Mr 

Lrprur 


(20) 


which  is  known  as  the  Reynolds  number.  Thus  the  normalized  form  of  Eqs.  (1-4) 
remains  invariant,  except  that  the  shear  stress  terms  have  to  be  multiplied  by 
the  factor  Re-1  (i.e.,  m  is  replaced  by  j_i*/Re).  As  an  example,  the  normalized 
form  of  txx  in  Eq.  (10)  becomes 


* 


T 


XX 


_  ★ 
2m 

3  Re 


av*  aw* 

.  *  .  *  ^ 

ay  az 


where  m*  is  given  from  the  Sutherland  law  as 


(21) 


* 

M 


(T*)3/2 


1  +  S/Tr 
T*  +  S/Tr 


(22) 


The  non-dimensional  heat  flux  term  in  the  energy  equation  takes  the  form 


— * 
q 


1  .  m  .  .  2>t  t  ai  t  aT  r 
(^r)( — r1  +  — T3  +  — 
ax  ay  az 


T-l  Re  '  Pr 


(23) 


Henceforth,  the  non-dimensional  equations  will  be  written  without  asterisks. 
Free  stream  conditions  are  used  in  this  report  as  the  reference  basis.  Hence, 
the  normalized  free  stream  values  of  density,  pressure,  and  temperature  will 
each  be  unity,  whereas  free  stream  velocity  components  will  be 


u_  =  Vtm  Cosa 

CD  00 


v  =  Vtm  Sina 

CD  CO 


(24) 


where  a  is  the  angle  of  attack  and  Mm  is  the  free  stream  Mach  number.  The 
normalized  free  stream  total  energy  per  unit  volume  from  Eq.  (5)  is 
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CD 


tm£ 

CD 
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+ 


1 

T-l 


(25) 
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A  new,  conservative,  finite  volume  scheme  is  developed  for  viscous 
(Navier-Stokes)  equations  as  an  extension  of  an  earlier  Euler  scheme  for 
unstructured  quadrilateral  grids.  Since  the  flux  balance  is  carried  out  on  an 
element  by  element  basis  (i.e.,  no  exterior  properties  are  utilized  for  a 
given  cell),  the  scheme  is  suitable  for  locally  embedded  unstructured  grids. 
The  viscous  and  inviscid  changes  for  a  given  time  step  can  be  decoupled  so 
that  cells  can  be  updated  independently  for  these  effects.  This  decoupling 
can  allow  substantial  savings  in  computer  resources  by  avoiding  computations 
of  viscous  terms  for  the  regions  where  they  are  negligible. 


3.3.1  Inviscid  Terms 


The  integration  scheme  for  inviscid  terms  has  been  previously  introduced 
in  Refs.  (9-12).  The  integration  scheme  is  based  upon  a  second  order 
accurate,  Lax-Wendrof f ,  finite  volume  scheme  due  to  Conservation 
variables  U  and  fluxes  F^ ,  are  all  stored  at  cell  vertices  and  flux  balance 
can  be  carried  out  for  each  cell  independently  for  a  given  time  step  At. 

Consider  the  flux  balance  for  a  cell  C  with  general  vertices  i,j,k,l,  as 
shown  in  Fig.  1,  which  results  in  the  distribution  formulas  for  the 
contribution  of  cell  change  to  its  corner  nodes 


sulc  •  4  t  alji  -  -  f-te  i 

sujc  -  i  [  »‘  *  ~SP  -  ^40  1 

s<4  -  ji™1*  raP  *  t45  1 

6u[c  -  \  t  aul  -  ^aP  ♦  1 


(26) 


The  superscript  i  signifies  the  inviscid  terms,  whereas  subscripts  like  iC 
signify  the  contribution  of  cell  C  to  node  i.  These  distribution  formulas 
allow  for  different  time  steps  At  and  cell  volumes  A  for  elements  adjoining  a 
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Figure  1.  An  Arbitrary  Cell  C  with  Nodes  i,j,k,l. 


common  node.  The  cell  change  AU1  in  the  above  relations  originates  from  a 
flux  balance  for  cell  C,  i.e. , 


^-AU1  =  F^(y  1-yi)  -  ^(x^Xi)  +  Fi(yk-yi)  -  G^x^x,) 


+  Fe<yj-yk>  -  Ge<xj~xk)  +  Fi<yi-yj)  -  Gs<xi-xj> 


(27) 


where  subscripts  s,  e,  n,  w  indicate  the  south,  east,  north,  west  fluxes, 
respectively.  Thus,  for  example 


F 


i 

s 


(28) 


The  flux  changes  AF  and  AS  are  due  to  second  order  changes  that  provide 
stability  to  the  algorithm  and  are  given  by 


AF  = 

AF^ 

<*! 

-  "g 

au 

)  AU 

AS  = 

AG1 

SF1 

(XS  iU 

-  y* 

au 

)  AU 
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The  flux  changes  involve  metrics  x^ ,  y^ ,  x^,  y^;  a  typical  example  of 
term  for  cell  C  is 


yn  -  yB 

Ayn 

yn  -  yB 


Ay 


yk  +  yi  -  yi  -  y-j 


ns  2 

yk  -*•  yx  -  yj  -  yi 


Ay 


ns 


The  Jacobian  matrix  Fy  is  given  by 


a  metric 


(30) 


(31) 
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-uv 
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u(F2i+u2-Ki| 

E^-+(1-Y)u2 

P 

uv(l-Y) 

uY 

The  Jacobian  matrix  Gy  is  similar  to  Fy  and  is  not  shown  here. 


3.3.2  Viscous  Terms 


A  new  scheme  is  presented  here  for  treating  the  viscous  terms  which,  like 
the  inviscid  scheme,  utilizes  only  the  properties  interior  to  a  given  cell  C. 
Only  the  first  order  contribution  from  the  viscous  terms  is  considered  here  to 
keep  the  computations  within  the  cell  c.  Note  that  the  only  reason  that  second 
order  terms  were  considered  for  the  inviscid  formulation  was  to  impart 
stability  to  the  numerical  algorithm,  since  a  first  order  Lax-Wendroff 
inviscid  scheme  is  inherently  unstable. 

The  flux  balance  around  the  boundary  dC  of  cell  C  of  Fig.  1  (See  Eqs. 
(26-27))  yields 


AU 

At 


?/,Fi 


dy  -  GLdx) 


2>C 


which  can  be  rewritten  as 


Z>C 


(33) 
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AU  =  AU1 


dy  -  Gvdx )  =  AU1  +  AUV 


(34) 


bC 

where  AU1  has  been  evaluated  in  the  previous  subsection  and  AUV  will  be 
evaluated  here.  This  equation  clearly  shows  that  the  viscous  and  inviscid 
flux  terms  can  be  decoupled.  It  is  apparent  from  Eq.  (26)  that  the  first 
order  contribution  to  change  at  ary  of  the  corner  nodes  is  one-fourth  of  the 
flux  balance  for  the  cell  C.  The  same  is  also  true  for  the  viscous  changes  and 
has  been  rigorously  proved  in  Ref.  (11)  for  transformed  coordinates.  Thus  the 
viscous  terms  can  be  handled  in  a  simple  manner  once  the  viscous  flux  balance 
for  the  cell  C  is  carried  out.  As  an  example  of  the  evaluation  of  AUV 
consider  the  contribution  from  x-momentum  equation.  Thus, 


AU(2)V  =  — 


At  _  bu 

jT{aw(yi  -  Vi)  <2^ 
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by 

bv 


bx  by 


&u  bv 

>s  -  Ps<xi  ■  xjH  —  +  —  >„> 


3  '  a y  &x 


(35) 


where  (2)  indicates  the  x-momentum  equation  contribution.  The  variable 
coefficients  a  and  p  are  given  by 


a  =  -lii- 

3  Re 


B  =  - 


Re 


(36) 


The  determination  of  these  coefficients  at  the  appropriate  faces  is  easy  to 
formulate,  e.g. , 


2  Re 


(ML  +  Mj) 


(37) 


where  p  is  evaluated  from  Eq.  (22).  The  determination  of  the  gradient  terms, 

I 

like  —  |w,  introduces  complications  if  one  insists  on  a  finite  volume 
approach.  One  has  to  either  employ  a  larger  staggered  stencil  (thereby  making 
the  scheme  unsuitable  for  unstructured  grids)  or  make  approximations  that 

the  metrics  and  areas  of  the  cells  surrounding  the  given  cell  are  nearly  equal 
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to  the  corresponding  quantities  in  cell  A  new  finite  volume  approach  is 

presented  here  for  unstructured  grids  in  which  such  approximations  are  avoided 

by  considering  auxiliary  staggered  cells  which  lie  completely  within  the  cell 

C.  These  interior  staggered  cells  are  shown  in  Fig.  2.  Thus  for  the 

determination  of  ^-|e  the  cell  jkns  is  used,  however  the  line  integration 

would  only  yield  the  gradient  at  the  centroid  ce  of  the  staggered  cell  rather 

than  at  the  edge  node  e.  However,  if  values  at  the  centroids  ce,  cn,  cw,  and 

cs  of  all  the  staggered  cells  can  be  determined,  then  one  can  extrapolate 

values  at  the  edge  nodes  e,  n,  w,  and  s.  As  an  example,  consider  the 
bu  iu 

evaluation  of  - —  and  - —  at  the  centroid  ce  by  applying  Green's  theorem  over 
2>x  by 

the  area  ijew.  The  line  integration  around  the  auxiliary  cell  yields^*^ 
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{  )  cs 


1  uw  +  ui 

— 2 — <vi  -  vw>  +  uc<yw  -  ye>  + 

ue  +  U1 

<  - 2 -  )(ye  ~  +  us(yj  "  yi)} 

1  . uw  +  ui 
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(38) 


(39) 


Here  the  values  u0,  u0,  etc.,  are  given  by  the  averages 

(40) 

(41) 


U8  "  7<Ui  + 

uc  =  i<ui  +  Uj  +  Uk  +  ul>  • 


Note  that  the  area  of  the  staggered  cells  has  been  approximated  here  as 
one-half  the  area  of  the  cell  C  under  consideration.  This  approximation  will 
not  hold  for  highly  skewed  cells.  However,  this  is  not  a  limitation  of  the 
overall  approach  since  the  actual  areas  of  the  staggered  cells  can  be  easily 
determined.  These  areas  will  have  to  be  determined  only  once  for  the  whole 
computat ion . 
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Figure  2.  Staggered  Cells  for  the  Computation  of  Edge 

o  U  | 

ijew  is  Used  for  — Cell  isnl  is  Used  for 
J>x  1  s 

A  biquadratic  extrapolation  can  be  used  to  determine  values  zs,  ?e,  zn, 
and  zw  once  some  z  values  at  the  centroids  of  the  auxiliary  cells  are 
determined.  These  extrapolations  yield 


Gradients,  Cell 

4U  I 

ixlw'  E 


zb  “  4(9zcs  +  zcn  "  3zce  “  3zcw^ 

zc  =  i<9zce  +  zcw  “  3 zc s  "  3zcn> 

zn  =  i<9zcn  +  zcs  -  3zce  "  3zcw> 

zw  =  ^<9zcw  +  zce  -  3zcs  -  3zcn>  * 


(42) 


The  computations  of  viscous  flux  changes  for  cell  C  can  then  be  determined  by 
using  an  equation  like  Eq.  (35).  The  viscous  contribution  to  changes  for 
nodes  i,j,k,l  is  then 


8U 


i»  j  f  k,  1 


(43) 
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3 . 4  Boundary  Conditions 


The  solution  of  the  fluid  dynamic  equations  is  determined  by  the  initial 
and  boundary  conditions.  The  initial  condition  for  the  cases  in  this  report 
is  typically  taken  as  that  corresponding  to  a  uniform  flow.  A  module  that 
utilizes  bilinear  interpolation,  based  upon  a  previous  solution  (on  a  somewhat 
different  grid),  has  also  been  developed  to  generate  initial  conditions  for  a 
given  new  grid  system.  Several  different  types  of  boundary  conditions,  such 
as  those  for  inflow,  outflow,  rigid  walls,  symmetry  walls,  etc.,  are  required 
for  computing  solutions.  Two  separate  inflow/outflow  boundary  conditions  are 
presented  here. 

3.4.1  Free-Slip  Rigid  Walls 

For  inviscid  flow,  the  appropriate  physical  condition  on  a  solid  surface 
is  that  there  be  no  flow  normal  to  the  surface,  or  equivalently  that  the  flow 
direction  be  tangential  to  the  wall.  The  application  of  the  free  slip  rigid 
wall  boundary  condition  is  explained  in  Fig.  3  for  a  node  i,  on  a  Bolid  wall 
which  makes  a  local  angle  a  with  the  x-axis. 


Figure  3.  Physical  Cells  and  their  Images  Adjacent  to  a  Wall. 
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B  after 


The  application  of  the  numerical  scheme  to  cells  A  and 
integrating  all  the  cells  in  the  domain  yields  the  following  initial  change  at 
node  i 

SUi  =  SUiA  +  SUiB  .  (44) 

The  predicted  change  at  node  i  is  taken  to  be  that  from  cells  A,B  and  their 
corresponding  mirror  images  A',B',  which  contribute  the  same  values,  i.e.. 


su?  = 

2SUi 

(45) 

The  superscript  p  indicates  predicted  change. 

If 

these  values 

are 

not 

corrected. 

then  the  wall  surface  would  be  a 

line 

of 

symmetry 

for 

all 

variables , 

including  of  course,  the  normal  component 

of 

the  velocity. 

The 

local  tangential  component  of  the  velocity  is  given  by 

Vfc  =  uCosa  +  vSina  ,  (46) 

and  only  this  is  used  to  reassign  new  velocity  components  along  the  coordinate 
directions,  i.e., 

u  =  VtCosa  ,  v  =  VtSina  .  (47) 

Thus,  the  corrected  values  for  the  changes  are 


S(pu)^  = 

(pV<_)*Cosa  - 

,  x  n 

(pu)  i 

Component  2 

(48) 

8(pv)?  = 

(pVt) *Sina  - 

(pv)i 

Component  3 

(49) 

SU? 

1 

suP 

Other  Components 

(50) 

where 

(pVt)*  =  [  (pu)?  +  S(pu)P]  Cosct  +  ((pv)J  +  S(pv)P)  Sinu  (51) 
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here  the  superscript  n  indicates  the  values  of  the  dependent  variables  at  the 
previous  time  level.  The  superscript  c  indicates  the  corrected  change  which, 
when  added  to  the  old  dependent  variables,  yields  the  values  at  the  new  time 
level  n+1. 


3.4.2  No-Slip  Rigid  Wall 

The  correct  constraint  at  a  viBcous  wall  is  the  no-slip  condition,  which 
implies  that  the  velocity  components  u,v  vanish  at  a  solid  boundary.  The 
temperature  can  either  be  prescribed  (possibly  isothermal  wall),  or  its 
gradient  can  be  set  to  specify  a  heat  flux  (possibly  adiabatic  wall).  The 
predicted  value  of  temperature  can  be  obtained  by  solving  the  caloric  equation 
of  state  and  accounting  for  the  predicted  change  in  density  and  energy  per 
unit  volume,  i.e., 

„  en+28e 

TP  =  (r-l)  - -  .  (52) 

pn+2Sp 

Note  that  the  corrected  value  of  density  is  taken  as 

pc  =  pn  +  25p  .  (53) 


The  corrected  value  of  temperature  at  the  boundary  node  is 
Tc  =  aT?  +  (l-a)Tw 


(54) 


where  a=l  for  an  adiabatic  wall  and  a=0  for  a  given  wall  temperature  Tw.  The 
corrected  value  of  energy  per  unit  volume  is  then 


€ 


c 


Tcpc 

T-l 


(55) 


3.4.3  Inflow/Outflow  Boundary  Conditions 

Two  separate  approaches  are  presented  here.  The  first  one  applies  the 
method  of  characteristics  for  Euler  equations  in  a  direction  normal  to  the 
local  grid  surfaces.  The  second  method  applies  the  usual  Riemann  invariants 
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for  isentropic  flows  at  these  boundaries.  The  former  approach  is  somewhat 
complicated  but  should  be  valid  for  cases  when  strong  shock  waves  intersect 
the  inflow/outflow  boundaries.  Both  the  approaches  neglect  the  viscous  terms 
at  the  exterior  flow  boundaries.  We  first  present  the  characteristic 

approach. 

Let  us  transform  the  flow  system  such  that  the  velocity  u^  corresponds  to 
that  entering  normal  to  the  grid  system  and  u^  corresponds  to  the  tangential 
direction.  This  is  shown  in  Fig  4.  The  following  relations  for  the 

transformation  of  velocities  hold 

u^  =  uSina  -  vCosa  (56) 

u^  =  uCosa  +  vSina  (57) 

with  the  inverse  relations 


u  =  u^Sina  +  u^Cosa  (58) 

v  =  u^Sina  -  u^Cosa  (59) 

where  a  is  the  angle  which  a  grid  boundary  surface  makes  with  x-axis. 

The  governing  equations  for  the  transformed  coordinates  are 


au  (  aF  t  i>G 
at  a^  an 


(60) 


where  only  the  inviscid  flux  vectors  are  considered.  It  is  assumed  that  the 
boundaries  are  far  enough  away  from  the  body  that  the  viscous  effects  are 
small.  Since  the  governing  equations  are  quasi-linear ,  they  can  be  written  as 


MJ 

it 


au 

as 


+  Gy 


au 

an 


o 


(61) 


The  Jacobian's  Fy  and  Gy  have  been  presented  in  Subsection  3.3.  We  will  now 
neglect  the  variations  along  the  tangential  direction,  i.e.,  set  au/an=0. 
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Figure  4.  Coordinate  Systems  for  Inflow/Outflow  Boundaries 


The  eigenvector  matrix  is  given  by 


V2 


+ 


V2 

2 


au 

Y-l 


au 

Y-l 


L  = 


-(u  + 


a 

Y-l 


) 


a 

Y-l 


u 


0 


1 


1 

1 

0 

1 


in  which  a  is  the  speed  of  sound  which  is  given  by 


,2  _  l£ 


The  corresponding  diagonal  matrix  for  the  eigenvector  matrix  L  is 


A  = 


u-a 

0 

0 

0 


0 

u+a 

0 

0 


0 

0 

u 

0 


Equation  (63)  represents  a  decoupled  system  of  equations  with  the 
change  in  characteristic  variables  given  by 

dQ  =  LdU 


or 


(65) 


(66) 


(67) 


differential 


(68) 
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(69) 


The  third  characteristic  variable  (third  row)  is  the  easiest  and  it  yields 

J  C 

pduT1=0  or  that  is  a  constant  along  the  characteristic  direction  —  =  u^. 
Other  characteristic  variables  yield  (after  some  simplification) 


dqj  = 

-  pa  duP  -  0 

Y-l  Y-l  S 

along  u^-  -  a 

(70) 

dq2  = 

aP  d  +  =  0 

T— 1  s  T— 1 

along  u^  +  a 

(71) 

dq3  = 

pdu^  =  0 

along  u^ 

(72) 

dp  a 

dq4  =  dp  =  0  along  u^  (73) 

Here  qj  (j=l,2,3,4)  are  regarded  to  be  the  elements  of  the  characteristic 
vector  variable  Q.  Assuming  the  coefficients  of  the  above  equations  to  be 
locally  frozen  (values  obtained  prior  to  updating),  these  equations  can  be 
easily  integrated.  The  frozen  values  are  shown  here  to  be  subscripted  by 
letter  c  and  the  characteristic  constants  are 
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<Jl 

=  p  - 

PcacuS 

(74) 

q2 

=  p  + 

Pcacui 

(75) 

^3 

=  un 

(76) 

<J4 

=  p  - 

2 

acp  . 

(77) 

Since  the  process  of  updating  (Un+^  =  Un  +SU)  involves  the  contributions 
to  nodes  from  the  interior  domain,  the  values  based  upon  Un+1  should  be  used 
for  characteristics  inferred  from  the  interior  computational  domain. 
Similarly,  the  characteristics  based  upon  external  conditions  should  either 
utilize  free  stream  values  or  values  based  upon  Un.  Wherever  these 
characteristic  variables  are  inferred  from,  the  primitive  variables  can  be 
computed  from  the  following  relations 


P 

qi+q2 

2 

(78) 

un 

=  q3 

(79) 

32-<Ji 

2Pcac 

(80) 

p 

qi+q2-2q4 

(81) 

For  each 

2a£ 

type  of  characteristic 

boundary  condition,  pc 

and  ac  are  determined 

from  Un. 

The  specific  boundary 

condition  prescription 

depends  upon  whether 

the  flow  field  is  subsonic  or  supersonic. 

For  supersonic  inflow  all  the  characteristic  information  comes  from  the 
exterior,  and  it  is  implemented  by  setting  the  changes  SU^  at  such  nodes  equal 
to  zero.  For  supersonic  outflow  all  the  characteristics  are  inferred  from  the 
interior,  and  the  changes  predicted  by  the  numerical  scheme  itself  can  be 
used.  Thus,  no  special  treatment  is  needed  for  supersonic  exit  boundaries. 
For  subsonic  inflow,  is  inferred  from  the  updated  Un  +  1  values  and  q2,  q3, 
are  based  upon  free  stream  values.  For  subsonic  exit  q3  is  computed  from 
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free  stream  values  (or  other  specified  exit  conditions),  whereas  q2,  q3,  q4 

are  computed  from  updated  Un+1  values. 

We  now  describe  the  alternate  inflow/outflow  boundary  conditions  based 
upon  Riemann  invariants.  For  the  subsonic  exit  the  outflow  pressure  p^  is 
specified.  For  supersonic  flow  with  some  imbedded  subsonic  regions  this 
pressure  is  inferred  from  the  supersonic  portion  of  the  flow  (at  the  place 
where  the  Mach  number  just  becomes  supersonic).  Let  p2,  u2,  v2,  p2  be  the 

flow  properties  at  some  node  interior  to  the  computational  domain.  The 
density  at  the  outflow  node  can  be  determined  from  the  constant  entropy 
condition 


Pi 


Pi  1/T 
P2  <— > 

P2 


(82) 


The  velocity  u-^  can  be  determined  from  the  R+  invariant  which  is  a  constant 
along  the  characteristic  =  u+a,  i.e.. 


U1  =  R+  “  7TT  ^rPl/pl 


(83) 


where 


R+  =  u2  +  A- VtP2/p2  .  (84) 

The  velocity  v^  is  determined  from  the  tangency  condition 
U1 

V1  =  v2<~ >  *  (85) 

Once  pj,  u^,  v^,  pj  are  known,  the  energy  can  be  determined  from  the 

caloric  equation  of  state  (5). 

For  the  subsonic  inflow  boundary,  total  pressure,  total  temperature  and 
v 

flow  angle  —  can  be  specified  while  the  Riemann  invariant  R_  can  be 
extrapolated  from  the  interior  of  the  flow.  Note  that 

R_  =  u  -  ;Ay VYp/p  (86) 

clx 

is  a  constant  along  the  characteristic  — —  =  u-a. 

dt 
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3 . 5  Stability  Analysis 


An  important  step  in  the  development  of  new  algorithms  is  the 
determination  of  time  step  restrictions  through  a  stability  analysis.  This 
section  explores  the  choice  of  a  time  step  for  the  explicit  schemes  presented 
in  Subsection  3.3  through  a  one-dimensional  convection-diffusion  model 
equation  by  utilizing  von  Neumann  stability  analysis.  Since  the  viscous  and 
inviscid  term  contributions  have  been  separated  out  in  Subsection  3.3, 
different  time  steps  can,  in  principle,  be  applied  to  these  terms.  The 
generalization  to  two-spatial  dimensions  is  also  carried  out  in  this  section. 

3.5.1  One-dimensional  Convection-dif fusion  Model 


The  one-dimensional  x-momentum  equation,  as  presented  in  Subsection  3.3, 
can  be  written  in  the  following  non-dimensional  form 


+ 


+  P) 


b  4m  <5U 

bx  3Re  bx  ^ 


(87) 


Considering  the  viscosity  p  to  be  a  constant  and  subtracting  off  the 
continuity  equation,  this  simplifies  to 


2>u  J>u  1  bp  4m  2>2u 

2>t  2>x  p  bx  3pRe  ^2 

The  model  equation  is  considered  to  be 
wave  speed  c,  i.e., 

J»U  2>U  J>2U 

—  +  c  —  =  v - 

at  ax  ax2 

where  one  can  regard  the  constant  v 
viscosity 


(88) 

similar  to  this  and  to  have  a  constant 

(89) 

to  be  proportional  to  the  kinematic 


v  = 

3pRe 

Discretizing  Eq.  (89)  according  to  a  Lax-Wendrof f ^ l*5*  scneme  yields 


(90) 
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u?+1  -  0? 

i.  1 


=  jUUi-i  -  Ui+1)  +  T(Ui_i  -  2UL  +  Ui.!)} 


+  EIU^  -  2UL  +  Ut+1) 


(91) 


where 


r  = 


cAt 

Ax2 


i9  the  CFL  (Courant-Friedrichs-Lewy )  number  and 


(92) 


Z  = 


vAt 


Ax" 


(93) 


is  the  reciprocal  of  the  cell  Reynolds  number,  which  shall  be  referred  to  here 
as  the  cell  diffusion  number.  Applying  the  standard  von  Neumann  stability 
analysis  to  Eq.  (91)  yields  the  following  constraint 

T2  +  2E  <  1  (94) 


This  constraint  is  a  combination  of  the  CFL  restriction 


T  <  1  (95) 

which  holds  for  inviscid  flows,  and  a  diffusion  stability  limitation 

z  \  (96> 

which  states  that  the  cell-Reynolds  number  has  to  be  greater  than  2  when  the 
convection  term  is  negligible.  On  account  of  Eq.  (95),  the  combined 
constraint  for  the  convection-diffusion  problem  can  be  revised  to  yield  a  more 
restrictive  form 


r  +  2Z  <  1 


(97) 


which  yields  the  following  expression  for  time  steps 


(17) 
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At  < 


TAX 

c+v /Ax£ 


(98) 


However,  since  the  viscous  and  inviscid  parts  cc.n  be  integrated  out  separately 
for  Ni  scheme,  there  is  no  need  to  implement  the  above  restrictive  form  on  the 
two  contributions.  The  inviscid  time  step  follows  from  the  CFL  restriction 
(Eq.  (95)) 


At,  < 


Ax 

u+a 


(99) 


where  u  is  the  flow  velocity  and  a  is  the  speed  of  sound,  whereas  the 
subscript  i  denotes  inviscid  time  step.  Note  that  the  constant  wave  speed  c 
has  been  replaced  by  the  maximum  eigenvalue,  u  +  a,  of  the  1-D  system.  The 
viscous  time  step  follows  from  the  diffusion  restriction 


At 


v 


< 


—  Ax2 

2  v 


3pReAx2E 

4m 


(100) 


where  the  subscript  v  indicates  the  viscous  time  step.  Hence  the  ratio  of  the 
two  time  steps  is  given  by 


Atv  _  3pReAx(u+a)E 

AtL  “  4ur 


(101) 


Utilizing  the  normalized  free  stream  values  for  these  quantities,  this  ratio 
of  the  two  stable  time  steps  can  be  shown  to  be 


Atv  3  r 

—  *  -  Re  Ax  Vr  ( 1+H(J0) 


(102) 


where  Y  is  the  ratio  of  specific  heats  and  is  the  free-stream  Mach  number. 
Thus  for  a  large  Reynolds  number  flow,  the  viscous  time  step  can  be  orders  of 
magnitude  larger  than  the  inviscid  time  step.  For  such  a  flow,  if  one 
utilizes  a  conventional  coupled  scheme  in  which  the  time  step  At  is  the 
minimum  of  At^  and  Atv  (or  a  hyperbolic  mean  as  suggested  by  Eq.  (98)),  then 
it  will  take  substantial  computer  resources  to  converge  on  the  viscous 
regions.  In  the  spirit  of  local  time  stepping,  the  viscous  and  inviscid 
contributions  can  be  updated  at  their  respective  time  restrictions  for  the 


29 


uncoupled  algorithm  developed  in  Subsection  3.3.  ThiB  will  allow  the 
algorithm  to  converge  significantly  faster  compared  to  conventional  uncoupled 
schemes . 

3.5.2  Time  Steps  in  2-D 

The  stability  limit  for  the  linearized  Euler  equations  from  the  analysis 
of  a  2-D  wave  equation  yields 

Ati  .  , _ 1 _  _ 1 _ 

TA  min  | uAy^-vAx^ | 4-aD^ '  | uAy^-vAx^ [ +aD^  ^  (1°3) 

where  A  is  the  area  of  a  cell  under  consideration  and 

=  Ax^  +  Ay^  .  (104) 

The  quantities  Ay^  =  Ayns,  Ay^  =  A yew,  etc.,  are  the  metrics  as  presented  in 
Subsection  3.3. 

The  generalization  for  the  viscous  time  steps  is  also  straight-forward, 

i.e.  , 


At 


v 


3pRe 

4u 


AE  min 


(105) 


3.5.3  Actual  Implementation  of  Time  Steps 

The  system  of  equations  becomes  numerically  stiff  for  very  high  Reynolds 
number.  It  is  clear  from  Eq.  (102)  that  for  such  a  case  the  viscous  and 
inviscid  time  steps  become  highly  disparate.  Since  these  time  step 
restrictions  have  been  derived  from  a  linearized  analysis,  it  is  no  surprise 
that  the  computations  become  unstable  for  non-linear  systems  of  equations  for 
high  Reynolds  number  flows.  The  stiffness  problem  is  addressed  in  this  study 
by  computing  a  maximum  time  step  for  each  cell 

Atm  =  F  min  ( At  £ ,  At v )  (106) 
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where  F  is  a  large  factor  that  keeps  At^  and  Atv  to  within  the  same  order  of 

magnitude;  a  maximum  value  of  F=10  will  accomplish  this.  The  inviscid  and 

viscous  time  steps  are  then  adjusted  as 

At^  =  max(At^,  Atm)  (107) 

Atv  =  max ( Atv,  Atm)  .  (108) 

Note  that  the  above  ad-hoc  recipe  for  time  steps  still  uses  different  values 
for  viscous  and  inviscid  terms,  but  limits  their  disparity  to  within  an  order 
of  magnitude  variation. 

3 . 6  Interactive  Block  Grid  Generator 


The  generation  of  initial  grids  for  complex  geometries  can  be  a  difficult 
task.  The  grid  generation  for  even  simple  flow  fields  with  multiple  embedded 
solid  objects  can  be  troublesome.  The  issue  of  the  complexity  of  the  grid 
generation  can  be  partially  simplified  by  considering  a  multiblock  gridding 
technique.  Subdomain  adaptation  or  automatic  division  of  cells  in  the  regions 
of  enhanced  gradients  of  selected  flow  field  variables,  can  be  used  to  further 
simplify  this  process. 

For  an  initial  grid  generation  an  interactive  block  grid  generator  (IBGG) 
has  been  developed  as  part  of  the  Phase  I  effort.  The  IBGG  subdivides  the 
flow  domain  into  blocks  which  are  simply  connected  regions.  Each  of  these 
blocks  is  associated  with  a  given  numerical  scheme.  These  blocks  are 
classified  as  Euler,  NS,  TLNS,  or  Interaction,  depending  upon  the  scheme  being 
utilized.  In  the  interaction  blocks  a  decision  is  locally  made  whether  to  use 
Euler  or  NS  equations. 

In  the  IBGG  approach  each  block  of  the  overall  configuration  initially 
has  an  independent  grid  structure.  The  topology  of  one  block  has  no  bearing 
on  the  rest  of  the  blocks,  except  that  there  must  be  a  node-to-node  matching 
across  the  interfaces  of  contiguous  blocks.  The  user  specifies  the  nodes  of 
the  boundaries  of  the  blocks  and  a  check  is  made  for  node-to-node  matching  at 
inter-block  interfaces.  Although  overlapping  interzonal  interfaces  provide 
maximum  geometrical  flexibility  for  complex  configurations,  such  overlapping 
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meshes  were  not  considered  for  this  research.  Overlapping  algorithms  are 
non-conservative  in  the  sense  that  interpolation  schemes  are  needed  to 
accomplish  intergrid  property  transfer.  Furthermore,  the  logic  of  coupling  an 
overlapping  grid  scheme  with  a  physics  and  subdomain  adaptation  would  have 
become  unduly  complicated. 

Figure  5  shows  the  first  step  in  the  interactive  block  grid  generation. 
The  geometry  corresponds  to  an  8%  circular  arc  cascade  with  a  flap  placed  near 
the  trailing  edge  of  the  arc.  The  user  has  selected  14  blocks.  The  circles 
correspond  to  the  user  selected  node  points  of  the  block  interfaces.  The  user 
has  clustered  the  node  points  near  the  solid  surfaces.  Only  part  of  the 
overall  block  system  is  shown.  As  shown  in  Fig.  6,  the  user  also  specifies 
the  numerical  schemes  to  be  used  on  individual  blocks.  Here  E  stands  for 
Euler  scheme,  V  for  a  viscous  scheme,  and  I  for  an  interaction  block. 

For  the  developed  1BGG ,  the  geometry  of  each  block  face  can  be  specified 
in  terms  of  cubic  polynomials.  Utilizing  the  interblock  boundary  point 
information,  the  interior  grid  of  each  of  the  blocks  is  generated  via  an 
algebraic  grid  generator.  This  step  also  determines  the  connectivity  arrays 
for  each  cell,  node  and  boundary  point  in  each  block.  Note  that  additional 
nodes  will  exist  at  the  time  of  assembly  of  the  overall  grid  when  the  points 
on  the  contiguous  boundaries  coincide.  These  multiply  defined  nodes  are 
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Figure  5.  User  Selected  Node  Points  for  14  Blocks  in  a  Computational 
Domain  of  an  8%  Circular  Arc  Cascade  with  a  Flap. 
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Figure  6.  User  Selected  Numerical  Schemes  for  the  14  Blocks  in  the 
Computational  Domain  of  an  8%  Circular  Arc  Cascade  with  a 
Flap. 


marked  for  deletion,  and  the  connectivity  arrays  of  the  interface  nodes  are 
examined  and  adjusted  for  consistency  in  the  data  structure.  Note  that  the 
total  number  of  cells  in  the  overall  grid  remains  invariant  as  the  blocks  are 
assembled.  Figure  7  shows  the  crude  base  mesh  generated  by  IBGG  after  the 
final  assembly  of  inter-block  interfaces.  At  this  point,  the  data  base 
forgets  about  the  blocks  and  all  the  alignments  are  in  terms  of  individual 
cells . 

After  the  generation  of  the  crude  assembly  all  the  cells  with  the  same 
numerical  algorithm  are  stored  together  in  the  computer  memory.  These  cells 
need  not  be  physically  contiguous.  At  this  stage  computations  can  be  carried 
out  on  the  crude  grid  since  all  the  connectivity  arrays  are  defined  and 
debugged.  Alternatively,  the  crude  assembly  can  be  filtered  through  an 
elliptic  grid  generator  to  move  the  nodes  appropriately.  This  filtering 
process  does  not  change  the  connectivity  arrays  but  rather  the  spatial 
coordinates  of  the  nodes.  An  elliptic  grid  generator  has  been  developed  that 
utilizes  the  pointer  system  of  the  unstructured  grid  connectivity  arrays.  In 
its  present  form,  it  fixes  the  locations  of  the  nodeB  on  the  physical 
boundaries,  but  allows  the  interior  nodes  (including  those  which  were 
originally  at  block  interfaces)  to  be  displaced.  Reference  (18)  presents  an 
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elliptic  grid  procedure  for  structured  multiblock  grids  where  special 
treatment  is  needed  at  the  block  interfaces.  No  such  special  treatment  is 
needed  for  the  overall  unstructured  data  base  utilized  here. 

Figure  8  shows  the  base  grid  after  the  IBGG  has  processed  the  previous 
crude  assembly  through  an  elliptic  grid  generator  which  had  no  forcing 
functions.  For  this  figure,  the  nodes  at  the  physical  boundaries  were  held 
fixed  and  the  absence  of  forcing  functions  caused  local  distortions  of  the 
grid  lines  near  these  boundaries.  Reference  (17)  presents  a  procedure  for  the 
determination  of  the  interior  grid  forcing  function  values  based  upon  their 
values  at  the  physical  boundaries.  Further  evaluation  of  forcing  functions 
will  be  carried  out  in  Phase  II. 

After  the  final  grid  assembly,  regions  of  cells  (irrespective  of  blocks) 
can  be  manually  subdivided  to  incorporate  pre-embedding.  Thus  cells  near 
solid  boundaries  can  be  manually  divided  to  capture  boundary  layers.  Such 
pre-embedding  can  result  in  enhanced  spatial  resolution  without  introducing 
grid  lines  in  the  regions  where  they  are  not  needed.  Each  cell  can  also  be 
spatially  divided  by  the  algorithm  to  incorporate  subdomain  adaptation. 
Physics  adaptation  may  be  carried  out  in  all  the  cells  which  previously 
belonged  to  interaction  blocks. 
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Figure  7.  Initial  Interior  Grid  Assembly  for  the  Computational  Domain 
of  an  8%  Circular  Arc  Cascade  with  a  Flap. 
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Figure  8.  Grid  Assembly  Produced  by  an  Elliptic  Grid  Generator  for  the 
Computational  Domain  of  an  8%  Circular  Arc  Cascade  with  a 
Flap. 

For  the  IBGG,  various  stages  of  the  grid  assembly  or  its  parts  can  be 
graphically  viewed  and  debugged  interactively.  The  user  could  request  the 
program  for  specific  changes.  For  example,  the  user  may  move  nodes  in  certain 
regions,  subdivide  or  fuse  meshes,  examine  or  query  specific  connectivities  of 
various  elements  (cells,  nodes,  boundary  points),  etc. 

3 . 7  Physics  Adaptation 

The  multiblock  SCRAP  technique  simplifies  physical  complexity  by 
associating  a  numerical  scheme  to  each  of  the  blocks.  Thus,  one  can  have 
viscous  blocks  where  NS  equations  are  solved  and  inviscid  blocks  where  one 
carries  out  the  solution  to  Euler  equations.  Such  a  scheme-splitting 
technique  can  efficiently  model  complex  physics  if  the  domains  of  viscous  and 
inviscid  regions  are  known  a  priori.  The  SCRAP  technique  addresses  this  issue 
by  considering  interaction  blocks  where  a  decision  is  locally  made  whether  to 
use  Euler  or  NS  equations.  This  dynamic  physics  switching  is  termed  here  as 
physics  (or  equation )  adaptation .  The  interaction  blocks  will  be  important 
for  the  cases  where  the  viscous/inviscid  interaction  region  develops  as  part 
of  the  numerical  solution. 
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There  are  numerous  advantages  for  utilizing  viscous/inviscid  interaction 
blocks  in  addition  to  separate  inviscid  and  viscous  blocks.  For  certain 
applications,  the  user  may  know  a  priori  exactly  where  the  inviscid  and 
viscous  parts  are  located.  For  example,  in  a  missile  configuration,  the 
"outer"  region  is  always  inviscid,  and  the  flow  near  the  surfaces  or  boundary 
layers  is  always  viscous.  For  such  distinct  regions,  computer  resources 
should  not  be  wasted  in  deciding  which  numerical  scheme  to  use.  Most  flow 
fields  also  possess  regions  where  it  is  not  clear  whether  to  use  a  viscous  or 
inviscid  formulation.  Such  regions  should  be  selected  as  interaction  blocks. 
Physics  adaptation  will  be  beneficial  for  these  intermediate  blocks  that  form 
the  buffer  between  viscous  and  inviscid  blocks.  If  such  interaction  blocks 
are  not  considered,  then  one  has  to  generally  treat  such  regions  as  being 
fully  viscous,  and  this  may  be  a  waste  of  resources.  In  addition,  if  viscous 
and  inviscid  blocks  are  directly  contiguous,  then  special  treatment  may  be 
need’d  at  block  interfaces  to  maintain  conservation.  It  can  be  shown  that 
such  treatment  is  not  needed  if  interaction  blocks  are  regarded  as  buffer 
zones  between  viscous  and  inviscid  blocks. 

As  shown  in  Subsection  3.3,  the  overall  algorithm  is  decoupled  into 
contributions  from  viscous  and  inviscid  terms.  This  means  that  Euler 
equations  are  solved  everywhere,  including  the  viscous  region.  This  is 
followed  by  the  determination  of  viscous  cells  in  the  interaction  blocks. 
These  cells  are  lumped  together  with  cells  in  the  viscous  regions.  Additional 
viscous  contributions  (corresponding  to  NS  or  TLNS  equations)  are  then 
computed  for  these  viscous  cells  and  added  to  the  inviscid  contributions.  The 
decoupled  nature  of  the  algorithm  allows  different  time  steps  to  be  taken  for 
viscous  and  inviscid  contributions.  This  is  in  line  with  the  concept  of  local 
time  stepping  for  carrying  out  solutions  to  steady  state  problems.  For  these 
steady  state  applications,  the  determination  of  viscous  cells  in  the 
interaction  blocks  can  be  carried  out  once  every  few  hundred  iterations.  This 
means  that  the  CPU  time  associated  with  the  process  of  selection  of  numerical 
schemes  will  be  small  compared  to  the  CPU  time  associated  with  the  integration 
of  cells.  For  unsteady  applications  where  the  viscous/inviscid  interaction 
regions  may  grow  significantly  as  part  of  the  numerical  solution,  physics 
adaptation  may  have  to  be  carried  out  more  frequently. 
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Physics  adaptation  is  incorporated  by  finding  the  maximum  shear  stress  in 
the  cells  of  viscous  and  interaction  blocks.  Those  cells  of  the  interaction 
block  for  which  the  shear  stress  is  more  than  1%  of  the  maximum  shear  stress 
are  regarded  as  viscous  cells.  Since  the  shear  stresses  only  change  by  a  very 
small  amount  near  the  interfaces  between  the  contiguous  viscous  and  inviscid 
cells,  a  special  treatment  is  not  needed  at  these  interfaces.  The  viscous  and 
inviscid  contributions  (respectively,  SuY  and  SU^  for  all  nodes  j)  are 
separated  out  for  each  node.  For  those  cells  which  are  only  regarded  as 
inviscid,  the  viscous  contributions  to  its  corresponding  nodes  is  regarded 
zero.  Obviously  at  the  interface  between  viscous  and  inviscid  cells,  some 
nodes  will  have  viscous  contributions  from  viscous  cells  and  the  same  nodes 
will  have  no  such  contribution  from  inviBcid  cells.  After  these  contributions 
are  summed,  the  process  of  updating  is  carried  out  by  adding  the  overall 
contributions  to  the  previous  values  of  the  dependent  variables. 

It  is  appropriate  to  point  out  that  the  placement  of  the  zonal  interface 
between  inviscid  and  viscous  solvers  is  very  important.  If  the  interface  is 
located  in  a  region  that  is  characterized  by  strong  flow  gradients  that  are 
not  modeled  by  Euler  equations,  then  there  would  appear  to  be  problems  with 
the  overall  solution.  This  is  not  really  a  shortcoming  of  the  zonal  gridding 
technique  but  rather  the  fact  that  important  physics  is  being  neglected  by 
using  the  inappropriate  mathematical  models.  Hence  care  must  be  exercised  in 
examining  the  overall  geometrical  configuration  for  the  initial  generation  of 
subdomain  grids  and  the  allocation  of  numerical  schemes  to  these  blocks.  It 
would  seem  logical  to  examine  the  relative  magnitude  of  the  stress  tensor  in 
the  inviscid  blocks,  after  a  specified  number  of  iterations,  to  confirm  the 
validity  of  the  inviscid  part.  If  the  stress  values  are  found  to  be 
substantial,  then  the  overall  algorithm  should  be  able  to  reassign  the 
numerical  schemes  in  corresponding  subdomains. 

As  part  of  a  feedback  from  ARO,  a  local  basis  for  the  selection  of 
viscous  cells  in  the  interaction  blocks  has  been  suggested.  Thus,  one  may 
compare  the  magnitude  of  viscous  terms  with  convection  (or  pressure)  terms 
locally  to  select  the  viscous  cells.  Although  the  global  approach  haB  worked 
fine  for  the  selected  model  problem,  its  disadvantage  is  that  two  passes  are 
needed  to  accomplish  the  objective.  In  the  first  pass,  maximum  viscous  stress 
is  determined  (from  viscous  and  interaction  blocks)  and  the  second  pass  tags 


37 


the  cells  in  the  interaction  region.  For  the  local  approach  only  a  single 
pass  o.ar  the  interaction  blocks  ceils  will  be  needed  to  tag  the  cells 
appropriately.  Since  the  suggested  approach  was  discussed  towards  the  end  of 
the  contract  period,  the  Principal  Investigator  (PI)  was  unable  to  evaluate 
it.  The  evaluation  and  utilization  of  a  local  approach  to  determine  viscous 
cells  will  be  carried  out  in  Phase  II. 

3 . 8  Subdomain  Adaptation 

The  SCRAP  technique  also  implements  subdomain  adaptation  in  which  the 
algorithm  automatically  implements  additional  spatial  resolution  by  locally 
dividing  cells  in  the  regions  of  enhanced  gradients  of  selected  flow  field 
variables.  The  details  of  this  process  have  appeared  in  Refs.  (9)— (12) .  This 
multi-variable  approach  defines  a  single  scalar  criterion  variable,  based  upon 
the  first  differences  of  multiple  components,  that  removes  the  effect  of 
inter-correlations  between  individual  components.  This  measure  allows  an 
unbiased  spread  of  data  for  the  cases  when  the  variabilities  in  different 
components  are  different  and  when  some  or  all  these  components  are  correlated. 
Since  the  single  scalar  variation  is  a  transformation  to  the  standardized 
form,  spatial  domains  characterized  by  different  kinds  of  variablities  can  be 
adapted  by  using  the  same  procedure.  Thus  the  same  approach  can  be  used  to 
adapt  in  the  viscous  and  inviscid  regions. 

The  cells  for  which  the  standardized  variation  exceeds  a  threshold  limit 
are  locally  subdivided.  The  cells  for  which  the  variations  remain  below 
another  smaller  threshold  value  are  allowed  to  merge.  First  differences  of 
density  were  used  to  adapt  in  the  inviscid  regions  near  shock  waves,  whereas 
the  first  differences  of  x-velocity  were  found  to  be  suitable  to  capture  the 
boundary  layer. 
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4.  PHASE  I  NUMERICAL  RESULTS 


This  section  describes  the  numerical  results  of  Phase  I.  The  potential  of 
the  SCRAP  technique  has  been  demonstrated  for  a  supersonic  flow  over  an  8% 
circular  arc  cascade.  This  model  problem  is  a  benchmark  case,  and  it  is  used 
to  test  the  current  approach.  Results  have  been  obtained  by  utilizing  various 
combinations  of  different  schemes  on  multiple  zones  (plus  varying  time  steps 
for  viscous  and  inviscid  schemes),  physics  adaptation  and  subdomain 
adaptation.  These  results  were  compared  with  those  utilizing  full  NS 
equations  on  a  single  fine  grid  and  yielded  essentially  the  same  answers  for 
the  same  boundary  conditions.  It  was  also  found  that  utilization  of  the 
characteristic  approach  versus  the  Riemann  invariant  approach  for  the  subsonic 
outflow  boundary  yielded  somewhat  different  results.  The  results  for  the 
subsonic  outflow  Riemann  invariant  approach,  utilizing  subdomain  and  physics 
adaptation,  compared  very  well  with  previous  computations  of  other 
researchers. 

The  variation  of  the  flow  properties  across  the  zonal  interfaces  was 
found  to  remain  independent  of  the  interface  itself.  The  fluid  properties 
also  remain  oblivious  to  those  interfaces  across  which  the  grid  changes 
geometrically  or  different  schemes  are  used.  The  technique  provides 
significant  savings  in  CPU  time  and  storage  over  conventional  structured  grid 
algorithms  employing  a  full  NS  solver.  An  increase  of  a  factor  of  about  50  in 
computational  speed  was  achieved  when  compared  to  a  conventional  full  NS 
solution  over  a  fine  grid  of  the  same  mesh  size  as  the  smallest  grid  size 
chosen  by  the  spatial  resolution  adaptation.  Further  savings  can  be  produced 
by  considering  additional  levels  of  subdomain  adaptation  to  capture  small 
scale  features.  This  demonstration  clearly  establishes  that  simultaneous 
physics  and  resolution  adaptation  works  for  the  sample  problem. 

4 . 1  Model  Problem 

The  potential  of  the  SCRAP  technique  can  be  demonstrated  by  considering  a 
model  problem  in  two  spatial  dimensions.  Figure  9  shows  the  geometry  and 
initial  fine  grid  for  such  a  model  problem.  The  distances  are  normalized  by 
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the  chord  length  for  this  8%  circular  arc  cascade.  The  free  stream  Mach 
number  is  1.4,  and  the  Reynolds  number,  based  upon  the  chord  length,  is 
23,000.  There  are  51  points  along  the  vertical  and  173  points  along  the 
horizontal  (8823  nodes).  This  model  problem  was  selected  because  the  solution 
for  both  Euler  and  NS  equations  is  well  known.  A  single  fine  grid  was 
initially  chosen  to  test  the  algorithm  and  so  that  comparisons  of  CPU  time  can 
be  made  with  other  adaptive  solutions.  Thus  there  are  initially  no  blocks  for 
the  solutions  of  Euler  and  full  NS  equations  in  this  single  grid. 

Supersonic  inlet  boundary  condition  is  applied  at  x=-l,  whereas 
character istic  outflow  condition  is  specified  at  x=2.0.  For  the  most  part, 
the  exit  condition  is  supersonic  with  a  subsonic  region  near  y=0.  The  top 
wall  is  considered  to  be  inviscid  solid  slip  wall  with  no  vertical  velocity. 
The  bottom  wall  is  a  slip  wall  for  Euler  equations,  and  a  no-slip  wall  for  NS 
equations.  Except  where  noted  otherwise,  the  boundary  condition  for  the 
subsonic  outflow  is  based  upon  the  characteristic  analysis  for  both  viscous 
and  inviscid  solutions. 


Figure  9. 


Geometry  and  Initial  Grid  (51X173)  For  Supersonic  Flow  Over 
an  8%  Circular  Arc  Cascade. 
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Figure  10.  Density  Contours  the  Euler  Solution  for  8%  Circular  Arc 

Cascade  on  a  51X173  Mesh  System.  The  Increment  for  Density 
Contours  is  0.05. 


4 . 2  Euler  Solution 


Figure  10  presents  the  Euler  solution  for  the  previous  geometry  and  grid. 
It  is  noted  that  the  scheme  has  good  shock  capturing  characteristics.  This 
solution  agrees  very  well  with  those  of  other  investigators. '  '  ’  The 
incident  shock  is  reflected  back  as  a  regular  reflected  shock  from  the  top 
boundary.  This  reflected  shock  intersects  with  the  trailing  edge  shock  and 
also  gets  reflected  at  the  lower  surface.  The  supersonic  outflow  boundary 
condition  is  clearly  appropriate  Bince  the  the  trailing  edge  shock  is  absorbed 
by  the  exit  boundary  rather  than  being  reflected.  The  total  CPU  time  for  this 
inviscid  solution  on  SSI's  Data  General  computer  was  72  minutes  when  the  RMS 
residuals  had  subsided  to  below  10~5. 

4 . 3  Full  Navier-Stokes  Solution 


Figure  11  represents  the  full  NS  solution  for  Re=23000  over  the  same  grid 
of  51x173  mesh  points.  The  NS  solution  is  distinctly  different  from  the  Euler 
solution.  The  boundary  layer  changes  the  flow  field,  and  this  means  that  the 
incoming  flow  now  encounters  an  effectively  larger  obstruction.  This  has  the 
effect  of  changing  the  flow  field  in  the  inviscid  part.  For  example  near  y=l, 
instead  of  a  regular  reflection,  a  Mach  disc  formation  is  observed.  The 
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reflected  shock  interacts  with  the  boundary  layer  and  results  in  shock 
wave/boundary  layer  interaction  which  causes  it  to  separate. 

Clearly,  additional  resolution  is  needed  to  resolve  the  separation 
region.  However,  for  the  sake  of  simple  CPU  time  comparisons  with  the  adapted 
cases  to  evaluate  the  proposed  scheme  in  Phase  I,  this  resolution  was  deemed 
adequate.  The  computed  solution  is  also  somewhat  different  from  that 
presented  in  Refs.  (15)  and  (21).  This  is  because  different  sorts  of  boundary 
conditions  were  used  for  the  subsonic  outflow  of  the  boundary  layer.  since 
substantial  time  and  effort  had  already  been  invested  in  this  fine  grid 
solution,  it  was  decided  to  retain  the  same  characteristic  subsonic  outflow 
boundary  condition  for  the  adapted  grid  solutions  presented  in  Subsections  4.4 
and  4.5.  A  solution  based  upon  Riemann  invariant  subsonic  outflow  boundary 
condition  was  only  carried  out  for  the  fully  adapted  solution,  and  it  is 
presented  in  Subsection  4.7.  The  latter  agrees  very  well  with  the  solutions 
in  the  cited  references. 

The  CPU  time  needed  to  converge  to  the  same  level  of  accuracy  as  the 
Euler  solution  of  Subsection  4.2  on  the  Data  General  machine  was  *ound  to  be 
178  minutes.  For  this  case  the  same  minimum  time  steps  were  used  for  the 
viscous  and  inviscid  contributions.  When  different  time  Bteps  were  used  for 
these  contributions,  the  total  convergence  time  was  165  minutes.  Thus  it  is 
found  that  using  different  time  steps  for  viscous  and  inviscid  contributions 
marginally  accelerates  the  convergence  process.  Since  the  solutions  for  these 
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Figure  11.  Density  Contours  for  the  Full  Navier-Stokes  Solution  for  8% 
Circular  Arc  Cascade  on  a  51X173  Mesh  System.  The  Increment 
for  Density  Contours  is  0.05. 
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two  cases  were  nearly  identical,  only  the  most  restrictive  time  step  solution 
is  presented  here.  The  CPU  time  for  the  full  NS  solution  is  about  a  factor  of 
two  longer  than  for  the  Euler  solution.  It  is  reasonable  to  assume  that  this 
ratio  is  Reynolds  number  dependent. 

4 . 4  Physics  Adaptation  Solution 

Figure  12  shows  the  grids  corresponding  to  the  three  blocks  that  were 
used  for  the  physics  adapted  solution.  The  composite  mesh  is  the  same  as  in 
the  previous  two  cases.  The  domain  is  split  into  Euler,  Interaction  and  NS 
blocks  and  only  part  of  the  Euler  block  is  shown.  The  viscous  block  is  chosen 
so  that  the  boundary  layer  will  be  captured.  The  interaction  block  will 
correspond  to  those  parts  of  the  boundary  layer  which  may  have  been  missed  in 
the  NS  block,  and  the  regions  where  the  reflected  shock  wave  interacts  with 
the  boundary  layer. 

Figure  13  shows  the  solution  corresponding  to  phyBics  adaptation 
procedure.  Those  cells  in  the  interaction  block  that  lie  near  the  region  of 
intersection  of  reflected  shock  wave  and  the  boundary  layer  are  automatically 
tagged  as  viscous  cells.  Except  for  some  minor  changes,  the  solution  is 
essentially  the  same  as  the  full  NS  solution. 


Figure  12.  Domain  Splitting  for  PhysicB  Adaptation  Solution  for  8% 
Circular  Arc  Cascade  on  a  Composite  Mesh  of  51X173  Points. 
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Figure  13.  Density  Contours  for  the  Physics  Adapted  Solution  for  8% 
Circular  Arc  Cascade  on  a  51X173  Initial  Grid.  Three 

Separate  Blocks  were  Used  for  NS,  Interaction  and  Euler 
Schemes.  The  Increment  for  Density  Contours  is  0.05. 


The  CPU  time  needed  to  converge  to  the  same  level  of  accuracy  as  the 
previous  two  solutions  was  found  to  be  82  minutes.  The  CPU  time  needed  for 
this  case  is  only  14%  more  than  that  required  for  the  full  Euler  solution. 
This  significant  saving  is  achieved  because  the  decision  basis  is  limited  only 
to  the  Interaction  block  and  because  different  time  steps  are  used  for  viscous 
and  inviscid  contributions. 


4 . 5  Physics  and  Subdomain  Adaptation  Solution 

Figure  14a  shows  the  spatially  adapted  grid  for  a  coupled  solution  of 
physics  adaptation  and  resolution  adaptation.  Two  levels  of  spatial  embedding 
are  carried  out  for  resolution  adaptation.  Note  that  the  base  grid  is  twice 
as  coarse  as  the  grid  used  for  previous  computations.  Furthermore,  the  finest 
mesh  is  half  as  fine  as  the  previous  case.  The  domain  splitting  for  the 
schemes  is  same  as  in  the  previous  computation. 

In  order  to  avoid  excessive  numerical  diffusion  of  the  boundary  layer  on 
this  initial  coarse  grid,  the  regions  in  the  boundary  layer  were  manually 
pre-embedded  to  a  single  level.  This  manual  subdivision  of  the  grid  is 
accomplished  as  part  of  the  solution  of  IBCG.  The  subdomain  adaptation  (both 
grid  division  and  fusion)  is  carried  out  automatically  after  a  specified 
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number  of  iterations,  and  as  is  apparent  fror..  tne  figure  it  can  be  achieved 
recursively.  There  are  11780  cells  and  9820  nodes  in  the  computational  domain 
for  the  final  grid.  First  differences  of  density  and  x-velocity  have  been 
used  to  accomplish  spatial  subdomain  adaptation. 

Figure  14b  shows  the  solution  corresponding  to  phy9ics/subdomain 
adaptation  procedure.  Additional  resolution  due  to  subdomain  adaptation 
allows  the  shock  and  boundary  layer  to  be  better  resolved  than  the  previous 
global  solutions.  The  solution  in  the  other  regions  is  very  nearly  the  same. 

The  CPU  time  needed  to  converge,  to  a  level  of  accuracy  of  RMS  errors 
being  less  than  10  5,  was  found  to  be  31  minutes.  This  means  that  the  fully 
adapted  solution  requires  only  17%  of  the  CPU  time  compared  to  that  of  the 
full  Navier-Stokes  solution.  However,  even  this  comparison  is  not  appropriate 
since  the  fully  adapted  solution  is  twice  as  fine  as  the  full  NS  solution. 
This  point  is  discussed  further  in  the  next  subsection. 
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Figure  14.  (a)  Final  Spatially  Adapted  Grid  Obtained  by  Utilizing  SCRAP 

Technique,  and  (b)  Density  Contours  on  this  Grid  with  an 
Increment  of  0.05. 
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4 . 6  CPU  Time  Comparisons 


Figure  15  summarizes  the  CPU  time  comparisons  for  various  cases.  The 
first  column  shows  the  actual  CPU  times  in  minutes  for  the  various  approaches. 
The  second  column  shows  the  normalized  computer  times,  where  the  Euler 
solution  is  regarded  as  the  reference  case.  The  third  column  indicates  how 
fast  a  particular  approach  is  in  reference  to  the  full  NS  solution.  The 
normalized  CPU  time  shows  that  full  NS  solution  is  2.5  times  slower  than  Euler 
solution,  Pereas  the  physics  adapted  solution  for  the  same  grid  takes  only 
14%  longer  than  the  Euler  equations  to  converge  for  the  model  problem. 

It  is  seen  that  the  fully  adapted  solution  converges  about  5.7  times 
faster  than  the  full  NS  solution.  However,  since  the  finest  resolution  of  the 
adapted  grids  is  half  as  much  as  the  global  grid,  a  rough  estimate  shows  that 
the  full  NS  solution  with  the  same  fine  grid  will  take  about  8  times  longer  to 
converge.  This  estimate  takes  into  account  that  there  are  four  times  as  many 
number  of  cells  with  twice  the  number  of  time  steps.  Furthermore,  physics 
adaptatioi  itself  is  about  2.5  times  faster.  This  translates  into  a  total 
additional  factor  of  20  w.r.t.  the  full  NS  solution.  Hence,  the  fully 

METHOD  ,  CPU  TIME  ,  NORMALIZED  ,  FACTOR 


(min) 
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72 
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PHYSICS  ADAPTATION 

82 

1.1 

2.2 

PHYSICS  S.  SUBDOMAIN 
ADAPTATION 

31 

0.4 

5.7* 

*  115  times  Foster  thon  corresponding  Full  NS  solution. 

Figure  15.  CPU  Time  Comparisons  for  Various  Approaches. 
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adapted  solution  is  about  two  orders  of  magnitude  faster  than  the 
corresponding  full  NS  solution.  With  additional  spatial  levels  of  adaptation, 
one  can  achieve  even  larger  factors.  For  3-D  computations,  this  saving  will 
be  even  more  substantial. 

4 . 7  Adapted  Solution  With  Alternate  Outflow  Boundary  Condition 

In  order  to  compare  the  solution  of  the  viscous  problem  with  those  in 
Refs.  (15)  and  (21),  the  fully  adapted  solution  was  carried  out  with  similar 
outflow  boundary  conditions.  The  subsonic  outflow  boundary  condition 
corresponds  to  the  extrapolation  of  Riemann  invariants  from  the  interior 
domain.  Pressure  was  inferred  from  the  supersonic  portion  of  the  boundary 
layer  and  was  held  fixed  for  the  subsonic  part.  Only  the  fully  adapted 
solution  was  repeated  because  it  runs  significantly  faster. 

The  initial  composite  grid  was  chosen  to  be  composed  of  25  x  25  points 
because  the  resolution  along  the  x-axiS  is  of  less  importance.  Furthermore, 
the  grid  in  the  boundary  layer  was  additionally  compressed  compared  to  the 
previous  computations.  The  scheme-splitted  blocks  were  similar  to  the 
previous  cases. 

Figure  16a  shows  the  final  adapted  grids  with  two  levels  of  spatial 
adaptation.  There  are  a  total  of  3852  nodes  and  4624  cells  in  this 
configuration. 

Figure  16b  shows  the  density  contours  for  the  grid  in  Fig.  16a.  It  is 
observed  that  the  boundary  layer  is  better  resolved  than  any  of  the  previous 
cases.  This  means  that  the  incomming  flow  observes  a  somewhat  smaller  bump, 
and  the  flow  is  less  choked.  The  Mach  disc  at  the  top  surface  is  less  severe, 
and  the  reflected  shock  wave  reaches  the  boundary  layer  aft  of  the  airfoil 
trailing  edge.  The  boundary  layer  re-attaches  itself  at  an  earlier  location, 
and  the  shock  reflection  emanating  from  the  boundary  layer  is  somewhat 
stronger. 

Figure  17  presents  the  contours  of  Mach  number.  These  agree  very  well 
with  the  computations  presented  in  Refs.  (15)  and  (21).  Figure  18  represents 
the  distribution  of  pressure  along  the  airfoil  surface;  the  agreement  with 
Refs.  (15)  and  (21)  is  again  very  good. 
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Figure  16.  (a)  Final  Spatially  Adapted  Grid  Obtained  by  Utilizing  SCRAP 

Technique,  and  (b)  Density  Contours  on  this  Grid  with  an 
Increment  of  0.05. 
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Figure  17.  Mach  Number  Contours  for  the  Fully  Adapted  Solution  for  an 
8%  Circular  Arc  Cascade.  The  Contour  Increments  are  0.05. 
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Figure  18.  Pressure  Distribution  Along  the  Lower  Channel  Wall  for  the 
8%  Circular  Arc  Cascade. 


The  results  presented  in  this  section  clearly  demonstrate  that 
simultaneous  physics  and  resolution  adaptation  works  for  the  sample  problem. 
The  technique  provides  significant  savings  in  CPU  time  and  storage  over 
conventional  structured  grid  algorithms  employing  a  full  NS  solver.  An 
increase  of  about  two  orders  of  magnitude  in  computational  speed  was  achieved 
when  compared  to  a  conventional  full  NS  solution  over  a  fine  grid  of  the  same 
mesh  size  as  the  smallest  grid  size  chosen  by  the  spatial  resolution 
adaptation. 
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5.  CONCLUSIONS  AND  SUMMARY 


The  Phase  I  effort  has  successfully  demonstrated  the  basic  feasibility  of 
simultaneous  physics  and  subdomain  adaptation  for  multiblock  grids.  The 
technique  provides  significant  savings  in  CPU  time  and  storage  over 
conventional  structured  grid  algorithms  employing  a  full  NS  solver  over  a  fine 
grid  of  the  same  mesh  size  as  the  smallest  grid  size  chosen  by  the  spatial 
resolution  adaptation.  The  SCRAP  technique  yields  essentially  the  same 
results  as  the  global  calculation  over  a  fine  grid  system.  The  results  were 
also  compared  with  previous  computations  and  were  found  to  be  in  good 
agreement.  This  successful  demonstration  of  proof  of  concept  in  Phase  I  has 
formed  the  basis  for  extension  to  3-D  in  Phase  II. 

The  results  presented  in  Section  4  clearly  demonstrate  that  simultaneous 
physics  and  resolution  adaptation  works  for  the  sample  problem.  An  increase 
of  about  two  orders  of  magnitude  in  computational  speed  has  been  achieved  when 
compared  to  a  conventional  full  NS  solution  over  fine  grids.  The  savings  in 
computer  storage  is  also  substantial  since  globally  fine  grids  are  not  needed. 
The  agreement  of  the  full'  \dapted  solution  with  globally  refined  solutions 
and  previous  computations  is  excellent.  The  features  that  enhance  the 
efficiency  of  the  technique  include  (1)  different  solvers  on  different  blocks 
and  their  dynamic  switching,  (2)  different  time  steps  for  viscous  and  inviscid 
contributions,  and  (3)  subdomain  adaptation.  Other  advantages  of  the  scheme 
include: 


•  ease  in  generating  grids  for  complex  geometries, 

•  ease  in  handling  complex  physics  of  specific  zones, 

•  accurate  modeling  of  features  like  shock  waves,  recirculation 
regions,  etc., 

•  greater  speed  compared  to  traditional  schemes, 

•  greater  accuracy  compared  to  traditional  schemes, 

•  comparatively  less  storage  required  for  same  accuracy. 
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The  advantages  for  the  SCRAP  technique  for  three  spatial  dimensions  will 
be  even  more  dramatic.  Formulation  of  the  data  management  connectivities  for 
3-D  will  be  carried  out  in  Phase  II.  This  work  will  produce  a  user-friendly 
software  package  that  will  alleviate  the  need  for  a  single  grid  and  allow 
concentration  of  computational  resources  in  modeling  complex  physics  only  in 
regions  of  importance.  This  technique  will  result  in  an  analysis  tool  that 
can  be  applied  in  carrying  out  flow  simulations  for  complex  and  irregular 
geometries  at  large  angles  of  attack.  Such  simulations  often  involve 
physically  complex  conditions  of  shock  waves,  boundary  layers,  separation  and 
recirculation  regions  and  their  interactions.  This  analysis  tool  will  more 
accurately  predict  the  resultant  forces  and  moments  on  practical  aerodynamic 
geometries.  This  will  lead  to  an  accurate  computation  of  the  trajectory  and 
design  of  guided  projectiles. 
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